// Copyright (c) Lawrence Livermore National Security, LLC and other VisIt
// Project developers.  See the top-level LICENSE file for dates and other
// details.  No copyright assignment is required to contribute to VisIt.

// ************************************************************************* //
//  File: avtLCSFilter.C
// ************************************************************************* //

#include <avtLCSFilter.h>
#include <avtLCSIC.h>

// For now use the avtLCSIC as the state does not need to be recorded
// for the FSLE. That is because currently the integration is being
// done step by step rather than in chunks. However, the code is set up
// to use avtIntegralCurveIC. Which if the integration is done in chucks
// will probably be more efficient.

//#include <avtIntegralCurveIC.h>
#define avtIntegralCurveIC avtLCSIC

#include <avtExtents.h>
#include <avtMatrix.h>
#include <avtParallel.h>
#include <avtCallback.h>

#include <avtOriginatingSource.h>
#include <vtkVisItScalarTree.h>

#include <VisItException.h>
#include <ImproperUseException.h>

#include <vtkMath.h>
#include <vtkUniformGrid.h>
#include <vtkRectilinearGrid.h>
#include <vtkStructuredGrid.h>
#include <vtkDoubleArray.h>
#include <vtkPointData.h>

#include <iostream>
#include <limits>
#include <cmath>

// ****************************************************************************
//  Method: avtLCSFilter constructor
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
//  Modifications:
//
//    Hank Childs, Wed Mar 28 08:36:34 PDT 2012
//    Initialize pathlines later.  Also set tolerances.
//
// ****************************************************************************

avtLCSFilter::avtLCSFilter() : seedVelocity(0,0,0)
{
    outVarRoot = std::string("operators/LCS/");

    replicateData = false;

    //outVarName ="operators/LCS/mesh";
    //doPathlines();
    //SetPathlines(atts.GetPathlines(),
    //              atts.GetPathlinesOverrideStartingTimeFlag(),
    //              atts.GetPathlinesOverrideStartingTime(),
    //              atts.GetPathlinesCMFE());
    //SetPathlines(true,false,0,PICS_CONN_CMFE);
    //SetPathlines(false,false,0,PICS_CONN_CMFE);

    // These initializations prevent harmless UMRs when we do our first
    // cache lookups.
    global_bounds[0] = global_bounds[2] = global_bounds[4] = 0;
    global_bounds[1] = global_bounds[3] = global_bounds[5] = 1;
    global_resolution[0] = global_resolution[1] = global_resolution[2] = 10;
    absTol = 1e-6;
    relTol = 1e-7;

    nDim = 3;
    auxIdx = LCSAttributes::None;
    nAuxPts = 1;
    auxSpacing = 0;

    numSteps = 0;
    fsle_dt = 0;
    fsle_ds = 0;

    cgTensor = LCSAttributes::Right;
    eigenComponent = LCSAttributes::Largest;

    minSizeValue =  std::numeric_limits<double>::max();
    maxSizeValue = -std::numeric_limits<double>::max();

    issueWarningForMaxStepsTermination = true;
    issueWarningForStepsize = true;
    issueWarningForStiffness = true;
    issueWarningForCriticalPoints = true;

    criticalPointThreshold = 1e-3;
}


// ****************************************************************************
//  Method: avtLCSFilter destructor
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
//  Modifications:
//
// ****************************************************************************

avtLCSFilter::~avtLCSFilter()
{
}


// ****************************************************************************
//  Method: avtLCSFilter::Create
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
// ****************************************************************************

avtFilter *avtLCSFilter::Create()
{
    return new avtLCSFilter();
}


// ****************************************************************************
//  Method: avtLCSFilter::SetAtts
//
//  Purpose:
//      Sets the state of the filter based on the attribute object.
//
//  Arguments:
//      a        The attributes to use.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
//  Modifications:
//
//    Hank Childs, Tue Mar 27 16:24:13 PDT 2012
//    Don't do pathlines if steady state is indicated.
//
//    Hank Childs, Wed Apr 11 11:35:16 PDT 2012
//    Add reverse flow.
//
// ****************************************************************************

void
avtLCSFilter::SetAtts(const AttributeGroup *a)
{
    atts = *(const LCSAttributes*)a;

    needsRecalculation =
      atts.ChangesRequireRecalculation(*(const LCSAttributes*)a);

    eigenComponent = atts.GetEigenComponent();
    clampLogValues = atts.GetClampLogValues();

    auxIdx = atts.GetAuxiliaryGrid();

    if( atts.GetOperationType() == LCSAttributes::IntegrationTime ||
        atts.GetOperationType() == LCSAttributes::ArcLength ||
        atts.GetOperationType() == LCSAttributes::AverageDistanceFromSeed )
    {
      nDim = 3;
      nAuxPts = 1;
      auxSpacing = 0;

      if( auxIdx != LCSAttributes::None )
      {
        avtCallback::IssueWarning("Requesting an auxiliary grid when none is needed. Ignoring request.");

        auxIdx = LCSAttributes::None;
      }
    }
    
    else if( auxIdx == LCSAttributes::None )
    {
      nDim = 3;

      nAuxPts = 1;
      auxSpacing = 0;
    }
    else if( auxIdx == LCSAttributes::TwoDim )
    {
      nDim = 2;
      nAuxPts = 4;
      auxSpacing = atts.GetAuxiliaryGridSpacing();
    }
    else if( auxIdx == LCSAttributes::ThreeDim )
    {
      nDim = 3;
      nAuxPts = 6;
      auxSpacing = atts.GetAuxiliaryGridSpacing();
    }

    cgTensor = atts.GetCauchyGreenTensor();
        
    int CMFEType = (atts.GetPathlinesCMFE() == LCSAttributes::CONN_CMFE
                    ? PICS_CONN_CMFE : PICS_POS_CMFE);

    SetPathlines(atts.GetPathlines(),
                 atts.GetPathlinesOverrideStartingTimeFlag(),
                 atts.GetPathlinesOverrideStartingTime(),
                 atts.GetPathlinesPeriod(),
                 CMFEType);

    SetIntegrationDirection(atts.GetIntegrationDirection());

    SetFieldType(atts.GetFieldType());
    SetFieldConstant(atts.GetFieldConstant());
    SetVelocitySource(atts.GetVelocitySource());

    SetIntegrationType(atts.GetIntegrationType());

    SetParallelizationAlgorithm(atts.GetParallelizationAlgorithmType(), 
                                atts.GetMaxProcessCount(),
                                atts.GetMaxDomainCacheSize(),
                                atts.GetWorkGroupSize());

    if (atts.GetIntegrationType() == LCSAttributes::DormandPrince)
    {
      if( atts.GetOperationType() ==  LCSAttributes::Lyapunov &&
            atts.GetTerminateBySize() )
        {
          EXCEPTION1(ImproperUseException,
                     "When performing FSLE the step size must be fixed. "
                     "Please select a different solver and "
                     "set the maximum time step.");
        }

        // For DoPri, the max time step is sent in to the PICS filter
        // as the max step length.
        double step;
        if (atts.GetLimitMaximumTimestep())
          step = atts.GetMaxTimeStep();
        else
          step = 0;

        SetMaxStepLength(step);
    }
    else
        SetMaxStepLength(atts.GetMaxStepLength());

    double absTol = 0.;
    bool doBBox = (atts.GetAbsTolSizeType() == LCSAttributes::FractionOfBBox);
    if (doBBox)
        absTol = atts.GetAbsTolBBox();
    else
        absTol = atts.GetAbsTolAbsolute();
    SetTolerances(atts.GetRelTol(), absTol, doBBox);

    SetTermination(atts.GetMaxSteps(),
                   atts.GetTerminateByDistance(),
                   atts.GetTermDistance(),
                   atts.GetTerminateByTime(),
                   atts.GetTermTime(),
                   atts.GetTerminateBySize(),
                   atts.GetTermSize());

    IssueWarningForAdvection(atts.GetIssueAdvectionWarnings());
    IssueWarningForBoundary(atts.GetIssueBoundaryWarnings());
    IssueWarningForMaxStepsTermination(atts.GetIssueTerminationWarnings());
    IssueWarningForStepsize(atts.GetIssueStepsizeWarnings());
    IssueWarningForStiffness(atts.GetIssueStiffnessWarnings());
    IssueWarningForCriticalPoints(atts.GetIssueCriticalPointsWarnings(),
                                  atts.GetCriticalPointThreshold());
}


// ****************************************************************************
//  Method: avtLCSFilter::Equivalent
//
//  Purpose:
//      Returns true if creating a new avtLCSFilter with the given
//      parameters would result in an equivalent avtLCSFilter.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
// ****************************************************************************

bool
avtLCSFilter::Equivalent(const AttributeGroup *a)
{
    return (atts == *(LCSAttributes*)a);
}


// ****************************************************************************
// Method: avtLCSFilter::SetVelocitySource
//
// Purpose: 
//   Sets the integral curve velocity source.
//
// Arguments:
//   vel : The velocity of the point.
//
//  Programmer: Allen Sanderson
//  Creation:   September 5, 2013
//
// ****************************************************************************

void
avtLCSFilter::SetVelocitySource(const double *p)
{
  seedVelocity.set(p);
}


// ****************************************************************************
//  Method: avtLCSFilter::GetInitialVelocities
//
//  Purpose:
//      Get the seed velocities out of the attributes.
//
//  Programmer: Allen Sanderson
//  Creation:   September 5, 2013
//
// ****************************************************************************

std::vector<avtVector>
avtLCSFilter::GetInitialVelocities(void)
{
    std::vector<avtVector> seedVels;

    seedVels.push_back( seedVelocity );

    return seedVels;
}


// ****************************************************************************
//  Method: avtLCSFilter::SetTermination
//
//  Purpose:
//      Sets the termination criteria for an integral curve.
//
//  Programmer: Allen Sanderson
//  Creation:   September 5, 2013
//
// ****************************************************************************

void
avtLCSFilter::SetTermination(int maxSteps_,
                              bool doDistance_, double maxDistance_,
                              bool doTime_,     double maxTime_,
                              bool doSize_,     double maxSize_)
{
    maxSteps = maxSteps_;
    doDistance = doDistance_;
    maxDistance = maxDistance_;
    doTime = doTime_;
    maxTime = maxTime_;
    doSize = doSize_;
    maxSize = maxSize_;
}


// ****************************************************************************
//  Method: avtLCSFilter::ModifyContract
//
//  Purpose:
//      Creates a contract the removes the operator-created-expression.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
// ****************************************************************************

avtContract_p
avtLCSFilter::ModifyContract(avtContract_p in_contract)
{
    avtDataRequest_p in_dr = in_contract->GetDataRequest();
    avtDataRequest_p out_dr = NULL;
    std::string var = in_dr->GetOriginalVariable();

    // This request is from an IC operator downstream.
    replicateData = in_contract->ReplicateSingleDomainOnAllProcessors();

//    in_contract->SetReplicateSingleDomainOnAllProcessors(true);
//    in_contract->SetOnDemandStreaming(false);
//    in_contract->GetDataRequest();
    in_dr->SetUsesAllDomains(true);

    if( strncmp(var.c_str(), "operators/LCS/", strlen("operators/LCS/")) == 0)
    {
        std::string justTheVar = var.substr(strlen("operators/LCS/"));

        outVarName = justTheVar;

        out_dr = new avtDataRequest(in_dr, justTheVar.c_str());
        //out_dr->SetDesiredGhostDataType(GHOST_NODE_DATA);
        //out_dr->SetDesiredGhostDataType(GHOST_ZONE_DATA);
    }
    else
    {
      outVarName = var;
      outVarRoot = "";
    }

    avtContract_p out_contract;

    if ( *out_dr )
        out_contract = new avtContract(in_contract, out_dr);
    else
        out_contract = new avtContract(in_contract);

    return avtPICSFilter::ModifyContract(out_contract);
}


// ****************************************************************************
//  Method: avtLCSFilter::UpdateDataObjectInfo
//
//  Purpose:
//      Tells output that we have a new variable.
//
//  Programmer: hchilds -- generated by xml2avt
//  Creation:   Mon Jan 10 07:15:51 PDT 2011
//
//  Modifications:
//    Brad Whitlock, Mon Apr  7 15:55:02 PDT 2014
//    Add filter metadata used in export.
//    Work partially supported by DOE Grant SC0007548.
//
// ****************************************************************************

void
avtLCSFilter::UpdateDataObjectInfo(void)
{
    avtDataAttributes  &in_dataatts =  GetInput()->GetInfo().GetAttributes();
    avtDataAttributes &out_dataatts = GetOutput()->GetInfo().GetAttributes();

    if( atts.GetSourceType() == LCSAttributes::RegularGrid )
    {
        out_dataatts.SetDataIsReplicatedOnAllProcessors(true);
        out_dataatts.SetDynamicDomainDecomposition(false);
    }

    //the outvarname has been assigned and will be added.
    //outVarName = "velocity";

    if (outVarRoot != "" && outVarName != "")
    {
        std::string fullVarName = outVarRoot + outVarName;

        out_dataatts.RemoveVariable(in_dataatts.GetVariableName());

        if (! out_dataatts.ValidVariable(fullVarName) )
        {
            //atts.AddVariable(outVarName.c_str());
            out_dataatts.AddVariable((fullVarName).c_str());
            out_dataatts.SetActiveVariable(fullVarName.c_str());
            //atts.SetTopologicalDimension(3);
            out_dataatts.SetVariableDimension(1);
            out_dataatts.SetVariableType(AVT_SCALAR_VAR);
            out_dataatts.SetCentering(AVT_NODECENT);
        }
    }
    else if (outVarName != "")
    {
      if( atts.GetOperationType() == LCSAttributes::EigenVector )
      {
        out_dataatts.SetVariableDimension(3);
        out_dataatts.SetVariableType(AVT_VECTOR_VAR);
      }
    }

    avtPICSFilter::UpdateDataObjectInfo();

    out_dataatts.AddFilterMetaData("LCS");
}


// ****************************************************************************
//  Method: avtLCSFilter::PostExecute
//
//  Purpose:
//      
//
//  Programmer: Hari Krishna
//  Creation:   December 5, 2011
//
//  Modifications:
//
// ****************************************************************************

void 
avtLCSFilter::PostExecute(void)
{
    avtPICSFilter::PostExecute();
}

// ****************************************************************************
//  Method: avtLCSFilter::PreExecute
//
//  Purpose:
//      Initialize data attributes for this filter and its base type (PICS).
//
//  Programmer: Hari Krishna
//  Creation:   December 5, 2011
//
//  Modifications:
//
//    Hank Childs, Jul  6 14:17:47 PDT 2012
//    Set resolution for Z to be 1 for 2D meshes.
//
// ****************************************************************************

void 
avtLCSFilter::PreExecute(void)
{
    SetActiveVariable(outVarName.c_str());
    GetSpatialExtents(global_bounds);

    // Cache the time and cycle at the start because for pathlines
    // they will change to the end.
    cycleCached = (int)   GetInput()->GetInfo().GetAttributes().GetCycle();
    timeCached  = (float) GetInput()->GetInfo().GetAttributes().GetTime();

    if (GetInput()->GetInfo().GetAttributes().GetSpatialDimension() == 2)
    {
        // we set them to 0->1 earlier and GetSpatialExtents only sets the
        // X and Y parts of the extents for 2D.
        global_bounds[4] = 0;
        global_bounds[5] = 0;
    }

    if(atts.GetUseDataSetStart() == LCSAttributes::Subset)
    {
        double* a = atts.GetStartPosition();
        global_bounds[0] = a[0];
        global_bounds[2] = a[1];
        global_bounds[4] = a[2];
    }

    if(atts.GetUseDataSetEnd() == LCSAttributes::Subset)
    {
        double* a = atts.GetEndPosition();
        global_bounds[1] = a[0];
        global_bounds[3] = a[1];
        global_bounds[5] = a[2];
    }

    const int* res = atts.GetResolution();
    global_resolution[0] = res[0];
    global_resolution[1] = res[1];
    global_resolution[2] = res[2];
    if (global_bounds[4] == global_bounds[5])
        global_resolution[2] = 1;

    double minResolution = std::numeric_limits<double>::max();

    double resX = 0., resY = 0., resZ = 0.;

    if (global_resolution[0] > 1)
    {
      double pcnt = 1.0 / (double) (global_resolution[0]-1);

      resX = (global_bounds[1] - global_bounds[0]) * pcnt; 

      minResolution = std::min( resX, minResolution );
    }
      
    if (global_resolution[1] > 1)
    {
      double pcnt = 1.0 / (double) (global_resolution[1]-1);

      resY = (global_bounds[3] - global_bounds[2]) * pcnt; 
 
      minResolution = std::min( resY, minResolution );
    }
      
    if (global_resolution[2] > 1)
    {
      double pcnt = 1.0 / (double) (global_resolution[2]-1);

      resZ = (global_bounds[5] - global_bounds[4]) * pcnt; 

      minResolution = std::min( resZ, minResolution );
    }

    if( atts.GetOperationType() == LCSAttributes::Lyapunov &&
        doSize && maxSize <= minResolution )
    {
        char str[1028];

        snprintf(str, 1028, "\nThe size limit for the FSLE is %f. "
                 "and is equal to or smaller than the resolution of the grid "
                 "(%f, %f, %f). ",
                 atts.GetTermSize(), resX, resY, resZ );

        avtCallback::IssueWarning(str);
        
//      EXCEPTION1(ImproperUseException, str );
    }

    avtPICSFilter::PreExecute();
}


// ****************************************************************************
//  Method: avtLCSFilter::Execute
//
//  Purpose:
//      Executes the LCS.  If we already have a cached version, then it
//      just returns that version.  If not, it calls PICS execute, which will
//      call our LCS set up routines via CreateIntegralCurveOutput.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

void
avtLCSFilter::Execute(void)
{
    avtDataTree_p dt = GetCachedDataSet();

    if (!needsRecalculation && *dt != NULL)
    {
        debug1 << "LCS: using cached version" << std::endl;

        SetOutputDataTree(dt);
    }
    else
    {
      debug1 << "LCS: no cached version, must re-execute" << std::endl;

      avtPICSFilter::Execute();

      std::vector<avtIntegralCurve *> ics;
      GetTerminatedIntegralCurves(ics);
      
      ReportWarnings( ics );
    }

#ifdef PARALLEL
    if (selectedAlgo == PICS_PARALLEL_OVER_DOMAINS && 
        atts.GetSourceType() == LCSAttributes::NativeMesh &&
        auxIdx == LCSAttributes::None )
      {
        char str[1028];

        snprintf(str, 1028,
                 "\n\nWhen utilizing a native mesh with mutliple domains "
                 "and no auxiliary grid it is not possible to compute "
                 "gradients across domain boundaries. For the best results "
                 "utilize an auxiliary grid or ghost cells." );

        avtCallback::IssueWarning(str);
      }
#endif


}

// ****************************************************************************
//  Method: avtLCSFilter::ContinueExecute
//
//  Purpose:
//      See if execution needs to continue.
//
//  Programmer: Allen Sanderson
//  Creation:   September 5, 2013
//
//  Modifications:
//
// ****************************************************************************

bool
avtLCSFilter::ContinueExecute()
{
    ++numSteps;

    if( atts.GetOperationType() == LCSAttributes::Lyapunov )
    {
      if( doSize )
      {
        std::vector<avtIntegralCurve *> ics;
        
        GetTerminatedIntegralCurves(ics);
        
        if (atts.GetSourceType() == LCSAttributes::NativeMesh)
        {
          if (NativeMeshIterativeCalc(ics) == true )
            return false;
          else
            return true;
        }
        else //if (atts.GetSourceType() == LCSAttributes::RegularGrid)
        {
          if (RectilinearGridIterativeCalc(ics) == true )
            return false;
          else
            return true;
        }
      }
    }

    return false;
}


// ****************************************************************************
//  Method: avtLCSFilter::GetInitialLocations
//
//  Purpose:
//      Tells the PICS filter where to place the initial seed locations.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

std::vector<avtVector>
avtLCSFilter::GetInitialLocations()
{
    seedPoints.clear();

    if (atts.GetSourceType() == LCSAttributes::NativeMesh)
    {
        GetInitialLocationsFromNativeMesh(GetInputDataTree());
    }
    else //if (atts.GetSourceType() == LCSAttributes::RegularGrid)
    {
        GetInitialLocationsFromRectilinearGrid();
    }

    return seedPoints;
}


// ****************************************************************************
//  Method: avtLCSFilter::GetInitialLocationsFromMesh
//
//  Purpose:
//      Walks through an AVT data tree and sets up the initial locations from
//      each point in the mesh.
//
//  Arguments:
//      inDT          A pointer to a data tree.  
 //
//  Programmer: Hank Childs
//  Creation:   December 5, 2011
//
// ****************************************************************************

void
avtLCSFilter::GetInitialLocationsFromNativeMesh(avtDataTree_p inDT)
{
    const double offset[3][6][3] = { { { 0., 0., 0.}, { 0., 0., 0.},
                                       { 0., 0., 0.}, { 0., 0., 0.},
                                       { 0., 0., 0.}, { 0., 0., 0.} },
                                     
                                     { {-1., 0., 0.}, { 1., 0., 0.},
                                       { 0.,-1., 0.}, { 0., 1., 0.},
                                       { 0., 0.,-1.}, { 0., 0., 1.} },
                                     
                                     { {-1., 0., 0.}, { 1., 0., 0.},
                                       { 0.,-1., 0.}, { 0., 1., 0.},
                                       { 0., 0.,-1.}, { 0., 0., 1.} } };
    if (*inDT == NULL)
        return;

    int nc = inDT->GetNChildren();

    if (nc < 0 && !inDT->HasData())
    {
        return;
    }

    if (nc == 0)
    {
        //
        // there is only one dataset to process
        //
        vtkDataSet *in_ds = inDT->GetDataRepresentation().GetDataVTK();
        size_t pts = in_ds->GetNumberOfPoints();

        size_t numberOfSeeds = seedPoints.size();
        size_t totalNumberOfSeeds = numberOfSeeds + pts*nAuxPts;
        seedPoints.resize( totalNumberOfSeeds );

        double base[3], point[3];
        
        // Insert the auxiliary points - the base point is not needed.
        for(size_t i = 0; i < pts; ++i)
        {
          in_ds->GetPoint(i, base);

          for(size_t a = 0; a < nAuxPts; ++a)
          {
            for(size_t l = 0; l < 3; ++l)
              point[l] = base[l] + auxSpacing * offset[auxIdx][a][l];

            seedPoints[numberOfSeeds++].set(point);
          }
        }
    }
    else
    {
        //
        // there is more than one input dataset to process
        // and we need an output datatree for each
        //
        for (int j = 0; j < nc; j++)
        {
            if (inDT->ChildIsPresent(j))
            {
                GetInitialLocationsFromNativeMesh(inDT->GetChild(j));
            }
        }
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::GetInitialLocationsFromRectilinearGrid
//
//  Purpose:
//      Created a series of seeds points based on the user specified mesh
//
//  Programmer: Hank Childs
//  Creation:   December 5, 2011
//
// ****************************************************************************

void
avtLCSFilter::GetInitialLocationsFromRectilinearGrid()
{
    const double offset[3][6][3] = { { { 0., 0., 0.}, { 0., 0., 0.},
                                       { 0., 0., 0.}, { 0., 0., 0.},
                                       { 0., 0., 0.}, { 0., 0., 0.} },
                                     
                                     { {-1., 0., 0.}, { 1., 0., 0.},
                                       { 0.,-1., 0.}, { 0., 1., 0.},
                                       { 0., 0.,-1.}, { 0., 0., 1.} },

                                     { {-1., 0., 0.}, { 1., 0., 0.},
                                       { 0.,-1., 0.}, { 0., 1., 0.},
                                       { 0., 0.,-1.}, { 0., 0., 1.} } };

    //compute total number of seeds that will be generated.
    size_t numberOfSeeds = 0;
    size_t totalNumberOfSeeds = nAuxPts *
      global_resolution[0] * global_resolution[1] * global_resolution[2];
    
    seedPoints.resize(totalNumberOfSeeds);

    double base[3], point[3];
        
    //add sample points by looping over in x,y,z
    for(int k = 0; k < global_resolution[2]; ++k)
    {
        double zpcnt = 0;

        if (global_resolution[2] > 1)
          zpcnt = ((double)k)/((double)global_resolution[2]-1);

        base[2] = global_bounds[4]*(1.0-zpcnt) + global_bounds[5]*zpcnt;
        
        for(int j = 0; j < global_resolution[1]; ++j)
        {
            double ypcnt = 0;

            if (global_resolution[1] > 1)
              ypcnt = ((double)j)/((double)global_resolution[1]-1);

            base[1] = global_bounds[2]*(1.0-ypcnt) + global_bounds[3]*ypcnt;
            
            for(int i = 0; i < global_resolution[0]; ++i)
            {
                double xpcnt = 0;

                if (global_resolution[0] > 1)
                  xpcnt = ((double)i)/((double)global_resolution[0]-1);

                base[0] = (global_bounds[0]*(1.0-xpcnt) +
                           global_bounds[1]*xpcnt);
                
                for(size_t a = 0; a < nAuxPts; ++a)
                {
                  for(size_t l = 0; l < 3; ++l)
                    point[l] = base[l] + auxSpacing * offset[auxIdx][a][l];

                  seedPoints[numberOfSeeds++].set(point);
                }
            }
        }
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::CreateIntegralCurve
//
//  Purpose:
//      Create an uninitialized integral curve.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

avtIntegralCurve*
avtLCSFilter::CreateIntegralCurve(void)
{
  if( doSize ) 
    return (new avtIntegralCurveIC());
  else
    return (new avtLCSIC());
}


// ****************************************************************************
//  Method: avtLCSFilter::CreateIntegralCurve
//
//  Purpose:
//      Create an integral curve with specific IDs and parameters.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

avtIntegralCurve*
avtLCSFilter::CreateIntegralCurve(const avtIVPSolver* model,
                                  avtIntegralCurve::Direction dir,
                                  const double& t_start,
                                  const avtVector &p_start,
                                  const avtVector& v_start, long ID)
{
    if (doPathlines)
    {
        if (dir == avtIntegralCurve::DIRECTION_BACKWARD)
            absMaxTime = seedTime0 - maxTime;
        else
            absMaxTime = seedTime0 + maxTime;
    }
    else
    {
        if (dir == avtIntegralCurve::DIRECTION_BACKWARD)
            absMaxTime = -maxTime;
        else
            absMaxTime = maxTime;
    }

    if( doSize )
    {
      // For now use the avtLCSIC as the state does not need to be
      // recorded for the FSLE. That is because currently the
      // integration is being done step by step rather than in
      // chunks. However, the code is set up to use
      // avtIntegralCurveIC. Which if the integration is done in chucks
      // will probably be more efficient.

      // unsigned char attr = avtStateRecorderIntegralCurve::SAMPLE_POSITION;
      // attr |= avtStateRecorderIntegralCurve::SAMPLE_TIME;
      // attr |= avtStateRecorderIntegralCurve::SAMPLE_ARCLENGTH;
      
      // return
      //   (new avtIntegralCurveIC(numSteps, doDistance, maxDistance, doTime, absMaxTime,
      //                        attr, model, dir, t_start, p_start, v_start, ID));
      return
        (new avtLCSIC(numSteps, doDistance, maxDistance, doTime, absMaxTime,
                      model, dir, t_start, p_start, v_start, ID));
    }
    else
    {
      return
        (new avtLCSIC(maxSteps, doDistance, maxDistance, doTime, absMaxTime,
                      model, dir, t_start, p_start, v_start, ID));
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::GetAllSeedsSentToAllProcs
//
//  Purpose:
//      
//
//  Programmer: Allen Sanderson
//  Creation:   August 5, 2015
//
// ****************************************************************************

bool
avtLCSFilter::GetAllSeedsSentToAllProcs(void)
{
  if (atts.GetSourceType() == LCSAttributes::NativeMesh)
  {
 #ifdef PARALLEL
    if (selectedAlgo == PICS_SERIAL)
      return true;
    else if (selectedAlgo == PICS_PARALLEL_OVER_DOMAINS)
      return false;
    // else if (selectedAlgo == PICS_PARALLEL_COMM_DOMAINS)
    //   return false;
    // else if (selectedAlgo == PICS_PARALLEL_MANAGER_WORKER)
    //   return false;
#else
    return true;
#endif
  }
  //  else if (atts.GetSourceType() == LCSAttributes::RegularGrid)
    return true;
}


// ****************************************************************************
//  Method: avtLCSFilter::CreateIntegralCurveOutput
//
//  Purpose:
//      Computes the LCS output (via sub-routines) after the PICS filter has
//      calculated the final particle positions.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

void 
avtLCSFilter::CreateIntegralCurveOutput(std::vector<avtIntegralCurve*> &ics)
{
  if( atts.GetOperationType() == LCSAttributes::Lyapunov && doSize )
  {
      if (atts.GetSourceType() == LCSAttributes::NativeMesh)
        CreateNativeMeshIterativeCalcOutput(ics);
      else //if (atts.GetSourceType() == LCSAttributes::RegularGrid)
        CreateRectilinearGridIterativeCalcOutput(ics);
  }
  else
  {
      if (atts.GetSourceType() == LCSAttributes::NativeMesh)
          NativeMeshSingleCalc(ics);
      else //if (atts.GetSourceType() == LCSAttributes::RegularGrid)
      {
          RectilinearGridSingleCalc(ics);

          // if( GetInputDataTree()->GetNChildren() )
          // {
          //   std::cerr << "Creating new data tree " << std::endl;
            
          //   avtDataTree_p outTree = MultiBlockDataTree( GetInputDataTree() );
            
          //   SetOutputDataTree(outTree);
          // }
      }
  }
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeLeftCauchyGreenTensor2D
//
//  Purpose:
//      Given a Gradient in place compute the Right Cauchy Green Tensor.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

void avtLCSFilter::ComputeLeftCauchyGreenTensor2D(double **j)
{
    double *j0 = j[0];
    double *j1 = j[1];
    
    // From the gradients compute the left Cauchy-Green Tensor J*J^T

    // j0[0]  j0[1] |  j0[0]  j1[0]
    // j1[0]  j1[1] |  j0[1]  j1[1]

    double a = j0[0]*j0[0] + j0[1]*j0[1];
    double b = j0[0]*j1[0] + j0[1]*j1[1];

    double d = j1[0]*j1[0] + j1[1]*j1[1];

    j0[0] = a;       j0[1] = b;
    j1[0] = b;       j1[1] = d;
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeLeftCauchyGreenTensor3D
//
//  Purpose:
//      Given a Gradient in place compute the Right Cauchy Green Tensor.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

void avtLCSFilter::ComputeLeftCauchyGreenTensor3D(double **j)
{
    double *j0 = j[0];
    double *j1 = j[1];
    double *j2 = j[2];
    
    // From the gradients compute the right Cauchy-Green Tensor J*J^T   

    // j0[0]  j0[1] j0[2] |  j0[0]  j1[0]  j2[0]
    // j1[0]  j1[1] j1[2] |  j0[1]  j1[1]  j2[1]
    // j2[0]  j2[1] j2[2] |  j0[2]  j1[2]  j2[2]

    double a = j0[0]*j0[0] + j0[1]*j0[1] + j0[2]*j0[2];
    double b = j0[0]*j1[0] + j0[1]*j1[1] + j0[2]*j1[2];
    double c = j0[0]*j2[0] + j0[1]*j2[1] + j0[2]*j2[2];

    double d = j1[0]*j1[0] + j1[1]*j1[1] + j1[2]*j1[2];
    double e = j1[0]*j2[0] + j1[1]*j2[1] + j1[2]*j2[2];
 
    double f = j2[0]*j2[0] + j2[1]*j2[1] + j2[2]*j2[2];

    j0[0] = a;       j0[1] = b;       j0[2] = c;
    j1[0] = b;       j1[1] = d;       j1[2] = e;
    j2[0] = c;       j2[1] = e;       j2[2] = f;
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeRightCauchyGreenTensor2D
//
//  Purpose:
//      Given a Gradient in place compute the Right Cauchy Green Tensor.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

void avtLCSFilter::ComputeRightCauchyGreenTensor2D(double **j)
{
    double *j0 = j[0];
    double *j1 = j[1];
    
    // From the gradients compute the right Cauchy-Green Tensor J^T*J

    // j0[0]  j1[0] | j0[0]  j0[1]
    // j0[1]  j1[1] | j1[0]  j1[1]

    double a = j0[0]*j0[0] + j1[0]*j1[0];
    double b = j0[0]*j0[1] + j1[0]*j1[1];

    double d = j0[1]*j0[1] + j1[1]*j1[1];

    j0[0] = a;       j0[1] = b;
    j1[0] = b;       j1[1] = d;
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeRightCauchyGreenTensor3D
//
//  Purpose:
//      Given a Gradient in place compute the Right Cauchy Green Tensor.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

void avtLCSFilter::ComputeRightCauchyGreenTensor3D(double **j)
{
    double *j0 = j[0];
    double *j1 = j[1];
    double *j2 = j[2];

    // From the gradients compute the right Cauchy-Green Tensor J^T*J   

    // j0[0]  j1[0]  j2[0] | j0[0]  j0[1] j0[2]
    // j0[1]  j1[1]  j2[1] | j1[0]  j1[1] j1[2]
    // j0[2]  j1[2]  j2[2] | j2[0]  j2[1] j2[2]


    double a = j0[0]*j0[0] + j1[0]*j1[0] + j2[0]*j2[0];
    double b = j0[0]*j0[1] + j1[0]*j1[1] + j2[0]*j2[1];
    double c = j0[0]*j0[2] + j1[0]*j1[2] + j2[0]*j2[2];

    double d = j0[1]*j0[1] + j1[1]*j1[1] + j2[1]*j2[1];
    double e = j0[1]*j0[2] + j1[1]*j1[2] + j2[1]*j2[2];
 
    double f = j0[2]*j0[2] + j1[2]*j1[2] + j2[2]*j2[2];

    j0[0] = a;       j0[1] = b;       j0[2] = c;
    j1[0] = b;       j1[1] = d;       j1[2] = e;
    j2[0] = c;       j2[1] = e;       j2[2] = f;
}


// ****************************************************************************
//  Method: avtLCSFilter::Jacobi2D
//
//  Purpose:
//      Computes the eigen vectors given a symmetric Jacobian.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

int avtLCSFilter::Jacobi2D(double **j, double *w)
{
    double *j0 = j[0];
    double *j1 = j[1];

    // Assume a symetric matrix
    // a b
    // b c

    double a = j0[0];
    double b = j0[1];
    double c = j1[1];

    double trace = (a + c) / 2.0;
    double det = a*c - b*b;
    double sqrtr = sqrt(trace * trace - det);

    // Order the largest first to match VTK
    w[0] = trace + sqrtr;
    w[1] = trace - sqrtr;

    return 1;
}

// ****************************************************************************
//  Method: avtLCSFilter::Jacobi2D
//
//  Purpose:
//      Computes the eigen vectors given a symmetric Jacobian.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

int avtLCSFilter::Jacobi2D(double **j, double *w, double **v)
{
    double *j0 = j[0];
    double *j1 = j[1];
    
    // Assume a symetric matrix
    // a b
    // b c

    double a = j0[0];
    double b = j0[1];
    double c = j1[1];

    double trace = (a + c) / 2.0;
    double det = a*c - b*b;
    double sqrtr = sqrt(trace * trace - det);

    // Order the largest first to match VTK
    w[0] = trace + sqrtr;
    w[1] = trace - sqrtr;

    // Use the largest value to get the first vector.
    avtVector vec0 = avtVector( -b, a-w[0], 0.0 );
    vec0.normalize();
    avtVector vec1 = avtVector( vec0.y, -vec0.x, 0.0 );

    v[0][0] = vec0.x;   v[0][1] = vec0.y; 
    v[1][0] = vec1.x;   v[1][1] = vec1.y; 


    // // Sorted in decreasing order.
    // double eigenvals[2];
    // double *eigenvecs[2];

    // double outrow0[2];
    // double outrow1[2];
    
    // eigenvecs[0] = outrow0;
    // eigenvecs[1] = outrow1;

    // vtkMath::JacobiN(j, 2, w, eigenvecs);

    // avtVector vec0 = avtVector( eigenvecs[0][0], eigenvecs[1][0], 0.0 );
    // vec0.normalize();
    // avtVector vec1 = avtVector( vec0.y, -vec0.x, 0.0 );

    // v[0][0] = vec0.x;   v[0][1] = vec0.y; 
    // v[1][0] = vec1.x;   v[1][1] = vec1.y; 

    return 1;
}

// ****************************************************************************
//  Method: avtLCSFilter::Jacobi3D
//
//  Purpose:
//      Computes the eigen vectors given a symmetric Jacobian.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

int avtLCSFilter::Jacobi3D(double **j, double *w)
{
    double *j0 = j[0];
    double *j1 = j[1];
    double *j2 = j[2];

    // Assume a symetric matrix
    // a b c
    // b d e
    // c e f

    double a = j0[0];
    double b = j0[1];
    double c = j0[2];
    double d = j1[1];
    double e = j1[2];
    double f = j2[2];

    // Now calculate the eigen values
    double x = ( a + d + f ) / 3.0f;  // trace

    a -= x;
    d -= x;
    f -= x;

    // Det / 2;
    double q = (a*d*f + b*e*c + c*b*e - c*d*c - e*e*a - f*b*b) / 2.0f;
    double r = (a*a + b*b + c*c + b*b + d*d + e*e + c*c + e*e + f*f) / 6.0f;

    double D = (r*r*r - q*q);
    double phi = 0.0f;

    if( D < std::numeric_limits<double>::epsilon())
      phi = 0.0f;
    else
    {
      phi = atanf( sqrt(D)/q ) / 3.0f;
      
      if( phi < 0 )
          phi += M_PI;
    }

    const double sqrt3 = sqrt(3.0f);
    const double sqrtr = sqrt(r);

    double sinphi = 0.0f, cosphi = 0.0f;
    sinphi = sinf(phi);
    cosphi = cosf(phi);

    // Sorted in decreasing order.
    double w0 = x + 2.0*sqrtr*cosphi;
    double w1 = x - sqrtr*(cosphi - sqrt3*sinphi);
    double w2 = x - sqrtr*(cosphi + sqrt3*sinphi);

    if( w0 >= w1 && w0 >= w2 )
    {
      w[0] = w0;

      if( w1 >= w2 ) { w[1] = w1;       w[2] = w2; }
      else           { w[1] = w2;       w[2] = w1; }
    }
    else if( w1 >= w0 && w1 >= w2 )
    {
      w[0] = w1;

      if( w0 >= w2 ) { w[1] = w0;       w[2] = w2; }
      else           { w[1] = w2;       w[2] = w0; }
    }
    else if( w2 >= w0 && w2 >= w1 )
    {
      w[0] = w2;

      if( w0 >= w1 ) { w[1] = w0;       w[2] = w1; }
      else           { w[1] = w1;       w[2] = w0; }
    }

    return 1;
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeEigenValues
//
//  Purpose:
//      Computes the eigen vectors given a Jacobian.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

void avtLCSFilter::ComputeEigenValues(vtkDataArray *jacobian[3],
                                      vtkDataArray *valArray)
{
    size_t nTuples = valArray->GetNumberOfTuples();

    if( auxIdx == LCSAttributes::TwoDim ||
        GetInput()->GetInfo().GetAttributes().GetSpatialDimension() == 2 )
    {
      for(size_t l = 0; l < nTuples; ++l)
      {
        double *input[2];
        input[0] = jacobian[0]->GetTuple3(l);
        input[1] = jacobian[1]->GetTuple3(l);

        if( cgTensor == LCSAttributes::Right )
          ComputeRightCauchyGreenTensor2D(input);
        else //if( cgTensor == LCSAttributes::Left )
          ComputeLeftCauchyGreenTensor2D(input);

        double eigenvals[2];
        Jacobi2D(input, eigenvals);

        if( eigenComponent == LCSAttributes::Largest )
          valArray->SetTuple1(l, eigenvals[0]);
        else // if( eigenComponent == LCSAttributes::Smallest )
          valArray->SetTuple1(l, eigenvals[1]);
      }
    }
    else
    {
      for(size_t l = 0; l < nTuples; ++l)
      {
        double *input[3];
        input[0] = jacobian[0]->GetTuple3(l);
        input[1] = jacobian[1]->GetTuple3(l);
        input[2] = jacobian[2]->GetTuple3(l);

        if( cgTensor == LCSAttributes::Right )
          ComputeRightCauchyGreenTensor3D(input);
        else //if( cgTensor == LCSAttributes::Left )
          ComputeLeftCauchyGreenTensor3D(input);

        double eigenvals[3];
        Jacobi3D(input, eigenvals);

        if( eigenComponent == LCSAttributes::Largest )
          valArray->SetTuple1(l, eigenvals[0]);
        else if( eigenComponent == LCSAttributes::Intermediate )
          valArray->SetTuple1(l, eigenvals[1]);
        else // if( eigenComponent == LCSAttributes::Smallest )
          valArray->SetTuple1(l, eigenvals[2]);
        }
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeEigenVectors
//
//  Purpose:
//      Computes the eigen vectors given a Jacobian.
//
//  Programmer: Allen Sanderson
//  Creation:   March 5, 2015
//
//  Modifications:
//
// ****************************************************************************

void avtLCSFilter::ComputeEigenVectors(vtkDataArray *jacobian[3],
                                       vtkDataArray *valArray,
                                       vtkDataArray *vecArray)
{
    double t0, t1, t2, sign, weight = atts.GetEigenWeight();

    if( eigenComponent == LCSAttributes::PosShearVector ||
        eigenComponent == LCSAttributes::PosLambdaShearVector )
      sign = +1;
    else
      sign = -1;
    
    size_t nTuples = valArray->GetNumberOfTuples();

    if( auxIdx == LCSAttributes::TwoDim ||
        GetInput()->GetInfo().GetAttributes().GetSpatialDimension() == 2 )
    {
      for(size_t l = 0; l < nTuples; ++l)
      {
        double *input[2];
        input[0] = jacobian[0]->GetTuple3(l);
        input[1] = jacobian[1]->GetTuple3(l);

        if( cgTensor == LCSAttributes::Right )
          ComputeRightCauchyGreenTensor2D(input);
        else //if( cgTensor == LCSAttributes::Left )
          ComputeLeftCauchyGreenTensor2D(input);

        // Sorted in decreasing order.
        double eigenval, eigenvals[2];
        double eigenvec[3], *eigenvecs[2];

        double outrow0[3];
        double outrow1[3];
        eigenvecs[0] = outrow0;
        eigenvecs[1] = outrow1;

        Jacobi2D(input, eigenvals, eigenvecs);

        if( eigenComponent == LCSAttributes::Largest )
        {
          eigenvecs[0][2] = 0;
          valArray->SetTuple1(l, eigenvals[0]);
          vecArray->SetTuple (l, eigenvecs[0]);
        }
        else if( eigenComponent == LCSAttributes::Smallest )
        {
          eigenvecs[1][2] = 0;
          valArray->SetTuple1(l, eigenvals[1]);
          vecArray->SetTuple (l, eigenvecs[1]);
        }
        else if( eigenComponent == LCSAttributes::PosShearVector ||
                 eigenComponent == LCSAttributes::NegShearVector )
        {
          eigenval = sqrt(eigenvals[0]) + sqrt(eigenvals[1]);

          t0 = sqrt( sqrt(eigenvals[0]) / eigenval );
          t1 = sqrt( sqrt(eigenvals[1]) / eigenval );

          // With the plus (minus) sign referring to the direction of
          // maximal positive (negative) shear in the frame of [ξ1,ξ0].
          eigenvec[0] = t0 * eigenvecs[1][0] + sign * t1 * eigenvecs[0][0];
          eigenvec[1] = t0 * eigenvecs[1][1] + sign * t1 * eigenvecs[0][1];
          eigenvec[2] = 0;

          valArray->SetTuple1(l, eigenval);
          vecArray->SetTuple (l, eigenvec);
        }
        // For calulating elliptic LCSs which are closed stationary
        // curves of the averaged strain. Solutions to this
        // variational problem are the closed orbits of one of
        // two parametrized vector-field families. 
        else if( eigenComponent == LCSAttributes::PosLambdaShearVector ||
                 eigenComponent == LCSAttributes::NegLambdaShearVector )
        {
          eigenval = weight * weight;
//        eigenval = eigenvals[1] * weight + eigenvals[0] * (1.0-weight);

          t0 = sqrt( (eigenvals[0]-eigenval) / (eigenvals[0]-eigenvals[1]) );
          t1 = sqrt( (eigenval-eigenvals[1]) / (eigenvals[0]-eigenvals[1]) );
          
          // With the plus (minus) sign referring to the direction of
          // maximal positive (negative) shear in the frame of [ξ1,ξ0].
          eigenvec[0] = t0 * eigenvecs[1][0] + sign * t1 * eigenvecs[0][0];
          eigenvec[1] = t0 * eigenvecs[1][1] + sign * t1 * eigenvecs[0][1];
          eigenvec[2] = 0;

          valArray->SetTuple1(l, eigenval);
          vecArray->SetTuple (l, eigenvec);
        }
      }
    }
    else //if( auxIdx == LCSAttributes::ThreeDim ||
         //    GetInput()->GetInfo().GetAttributes().GetSpatialDimension() == 3 )
    {
      for(size_t l = 0; l < nTuples; ++l)
      {
        double *input[3];
        input[0] = jacobian[0]->GetTuple3(l);
        input[1] = jacobian[1]->GetTuple3(l);
        input[2] = jacobian[2]->GetTuple3(l);

        if( cgTensor == LCSAttributes::Right )
          ComputeRightCauchyGreenTensor3D(input);
        else //if( cgTensor == LCSAttributes::Left )
          ComputeLeftCauchyGreenTensor3D(input);

        // Sorted in decreasing order.
        double eigenval,     eigenvals[3];
        double eigenvec[3], *eigenvecs[3];

        double outrow0[3];
        double outrow1[3];
        double outrow2[3];
        eigenvecs[0] = outrow0;
        eigenvecs[1] = outrow1;
        eigenvecs[2] = outrow2;

        vtkMath::Jacobi(input, eigenvals, eigenvecs);

        // vtkMath::Jacobi and vtkMath::JacobiN return the vectors in
        // a column ordering. That is the first vector is in:
        // eigenvecs[0][0] eigenvecs[1][0] eigenvecs[2][0] and that
        // output[0] output[1] output[1] does NOT contain the first
        // eigen value but the first component of each eigen vector.
            
        if( eigenComponent == LCSAttributes::Largest )
        {
          valArray->SetTuple1(l, eigenvals[0]);
          vecArray->SetTuple3(l, eigenvecs[0][0], eigenvecs[1][0], eigenvecs[2][0]);
        }
        else if( eigenComponent == LCSAttributes::Intermediate )
        {
          valArray->SetTuple1(l, eigenvals[1]);
          vecArray->SetTuple3(l, eigenvecs[0][1], eigenvecs[1][1], eigenvecs[2][1]);
        }
        else if( eigenComponent == LCSAttributes::Smallest )
        {
          valArray->SetTuple1(l, eigenvals[2]);
          vecArray->SetTuple3(l, eigenvecs[0][2], eigenvecs[1][2], eigenvecs[2][2]);
        }
        else if( eigenComponent == LCSAttributes::PosShearVector ||
                 eigenComponent == LCSAttributes::NegShearVector )
        {
          eigenval = sqrt(eigenvals[0]) + sqrt(eigenvals[2]);

          t0 = sqrt( sqrt(eigenvals[0]) / eigenval );
          t2 = sqrt( sqrt(eigenvals[2]) / eigenval );

          // With the plus (minus) sign referring to the direction of
          // maximal positive (negative) shear in the frame of [ξ2,ξ0].
          eigenvec[0] = t0 * eigenvecs[0][2] + sign * t2 * eigenvecs[0][0];
          eigenvec[1] = t0 * eigenvecs[1][2] + sign * t2 * eigenvecs[1][0];
          eigenvec[2] = t0 * eigenvecs[2][2] + sign * t2 * eigenvecs[2][0];

          valArray->SetTuple1(l, eigenval);
          vecArray->SetTuple (l, eigenvec);
        }
        else if( eigenComponent == LCSAttributes::PosLambdaShearVector ||
                 eigenComponent == LCSAttributes::NegLambdaShearVector )
        {
          eigenval = weight * weight;
//        eigenval = eigenvals[1] * weight + eigenvals[0] * (1.0-weight);

          t0 = sqrt( (eigenvals[0]-eigenval) / (eigenvals[0]-eigenvals[2]) );
          t2 = sqrt( (eigenval-eigenvals[2]) / (eigenvals[0]-eigenvals[2]) );
          
          // With the plus (minus) sign referring to the direction of
          // maximal positive (negative) shear in the frame of [ξ2,ξ0].
          eigenvec[0] = t0 * eigenvecs[0][2] + sign * t2 * eigenvecs[0][0];
          eigenvec[1] = t0 * eigenvecs[1][2] + sign * t2 * eigenvecs[1][0];
          eigenvec[2] = t0 * eigenvecs[2][2] + sign * t2 * eigenvecs[2][0];

          valArray->SetTuple1(l, eigenval);
          vecArray->SetTuple (l, eigenvec);
        }
      }
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::ComputeLyapunovExponent
//
//  Purpose:
//      Computes the FTLE or FLLE given a Jacobian. Which is the
//      following: log of the square root of the maximum eigen value
//      of the right Caughy-Green Tensor with the result divided by
//      the time or distance.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
//  Modifications:
//
//    Hank Childs, Fri Sep  7 15:47:12 PDT 2012
//    Convert calculation to double precision, which prevents a "cliff" from
//    too large epsilon associated with float.
//
// ****************************************************************************

void avtLCSFilter::ComputeLyapunovExponent(vtkDataArray *jacobian[3],
                                           vtkDataArray *expArray)
{
    bool takeLog = doTime || doDistance;

    size_t nTuples = expArray->GetNumberOfTuples();

    // The base value is used to clamp the log values to be only
    // positive or both positive and negative.
    double baseValue, denominator = 1.0;

    // Clamp only if taking the log.
    if (takeLog && clampLogValues == true )
      baseValue = 1.0;
    else
      baseValue = -std::numeric_limits<double>::max();

    if( doTime )
      denominator /= maxTime;
    else if( doDistance )
      denominator /= maxDistance;

    if( auxIdx == LCSAttributes::TwoDim ||
        GetInput()->GetInfo().GetAttributes().GetSpatialDimension() == 2 )
    {
      for(size_t l = 0; l < nTuples; ++l)
      {
        // if( (doTime &&
        //      (fabs(expArray->GetTuple1(l) - absMaxTime) > FLT_MIN)) ||
        //     (doDistance &&
        //      (expArray->GetTuple1(l) < maxDistance)) )
        //   expArray->SetTuple1(l, 0);
        // else
        {
          double *input[2];
          input[0] = jacobian[0]->GetTuple3(l);
          input[1] = jacobian[1]->GetTuple3(l);

          if( cgTensor == LCSAttributes::Right ) 
            ComputeRightCauchyGreenTensor2D(input);
          else //if( cgTensor == LCSAttributes::Left )
            ComputeLeftCauchyGreenTensor2D(input);

          // Get the eigen values.
          double  eigenvals[2];
          Jacobi2D( input, eigenvals );

          double lambda = baseValue;

          if( eigenComponent == LCSAttributes::Largest )
            lambda = sqrt( std::max( lambda, eigenvals[0] ) );
          else // if( eigenComponent == LCSAttributes::Smallest )
            lambda = sqrt( std::max( lambda, eigenvals[1] ) );

          if( takeLog )
            lambda = log( lambda );

          lambda *= denominator;

          expArray->SetTuple1(l, lambda);
        }
      }
    }
    else //if( auxIdx == LCSAttributes::ThreeDim ||
         //    GetInput()->GetInfo().GetAttributes().GetSpatialDimension() == 3 )
    {
      for(size_t l = 0; l < nTuples; ++l)
      {
        // if( (doTime &&
        //      (fabs(expArray->GetTuple1(l) - absMaxTime) > FLT_MIN)) ||
        //     (doDistance &&
        //      (expArray->GetTuple1(l) < maxDistance)) )
        //   expArray->SetTuple1(l, 0);
        // else
        {
          double *input[3];
          input[0] = jacobian[0]->GetTuple3(l);
          input[1] = jacobian[1]->GetTuple3(l);
          input[2] = jacobian[2]->GetTuple3(l);

          if( cgTensor == LCSAttributes::Right ) 
            ComputeRightCauchyGreenTensor3D(input);
          else //if( cgTensor == LCSAttributes::Left ) 
            ComputeLeftCauchyGreenTensor3D(input);

          // Get the eigen values.
          double eigenvals[3];
          Jacobi3D( input, eigenvals );

          double lambda = baseValue;

          if( eigenComponent == LCSAttributes::Largest )
            lambda = sqrt( std::max( lambda, eigenvals[0] ) );
          else if( eigenComponent == LCSAttributes::Intermediate )
            lambda = sqrt( std::max( lambda, eigenvals[1] ) );
          else // if( eigenComponent == LCSAttributes::Smallest )
            lambda = sqrt( std::max( lambda, eigenvals[2] ) );

          if( takeLog )
            lambda = log( lambda );

          lambda *= denominator;

          expArray->SetTuple1(l, lambda);
        }
      }
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::ReportWarnings() 
//
//  Purpose:
//      Reports any potential integration warnings
//
//  Programmer: Allen Sanderson
//  Creation:   20 August 2013
//
//  Modifications:
//
// ****************************************************************************

void
avtLCSFilter::ReportWarnings(std::vector<avtIntegralCurve *> &ics)
{
    int numICs = (int)ics.size();

    int numAdvection = 0;
    int numBoundary = 0;

    int numEarlyTerminators = 0;
    int numStepSize = 0;
    int numStiff = 0;
    int numCritPts = 0;

    if (DebugStream::Level5())
    {
        debug5 << "::ReportWarnings " << ics.size() << endl;
    }

    // Loop through all the IC for warnings.
    if (doSize)
    {
        for (int i = 0; i < numICs; ++i)
        {
            avtIntegralCurveIC *ic = dynamic_cast<avtIntegralCurveIC*>(ics[i]);

            bool badSize = (doSize && ic->GetTime() < maxSize);

            if (ic->CurrentVelocity().length() <= criticalPointThreshold)
              ++numCritPts;

            if (ic->TerminatedBecauseOfMaxSteps())
              ++numEarlyTerminators;
            
            if (ic->status.StepSizeUnderflow() && badSize)
              ++numStepSize;

            if (ic->EncounteredNumericalProblems())
              ++numStiff;

            if (ic->status.ExitedSpatialBoundary())
              ++numBoundary;

            if (badSize)
              ++numAdvection;         
        }
    }
    else // if (doTime || doDistance)
    {
        for (int i = 0; i < numICs; ++i)
        {
            avtLCSIC *ic = dynamic_cast<avtLCSIC*>(ics[i]);

            bool badTime = (doTime && (fabs(ic->GetTime() - absMaxTime) > FLT_MIN));
            bool badDistance = (doDistance && (ic->GetDistance() < maxDistance));

            if (ic->CurrentVelocity().length() <= criticalPointThreshold)
              ++numCritPts;

            if (ic->TerminatedBecauseOfMaxSteps())
              ++numEarlyTerminators;
            
            if (ic->status.StepSizeUnderflow() && (badTime || badDistance))
              ++numStepSize;

            if (ic->EncounteredNumericalProblems())
              ++numStiff;

            if (ic->status.ExitedSpatialBoundary())
              ++numBoundary;

            if (badTime || badDistance)
              ++numAdvection;
        }
    }

    char str[4096] = "";

    if (issueWarningForAdvection)
    {
        SumIntAcrossAllProcessors(numAdvection);

        if (numAdvection)
        {
          snprintf(str, 4096,
                   "%s\n%d of your integral curves terminated before they reached "
                   "the maximum advection criteria.  This may be indicative of your "
                   "time, distance, or size criteria being too large or the curve leaving the domain."
                   "  Note that this message does not mean that an error has occurred; it simply "
                   "means that VisIt stopped advecting particles before they reached the maximum.\n",
                   str, numAdvection);
        }
    }

    if ((doDistance || doTime || doSize) && issueWarningForBoundary)
    {
        SumIntAcrossAllProcessors(numBoundary);
        if (numBoundary > 0)
        {
            snprintf(str, 4096, 
                     "%s\n%d of your integral curves exited the spatial domain. This will all likelihood, "
                     "skew the results of any relative gradient based (Eigen value/vector) calculations.\n", str, numBoundary);
        }
    }

    if ((doDistance || doTime || doSize) && issueWarningForMaxStepsTermination)
    {
        SumIntAcrossAllProcessors(numEarlyTerminators);
        if (numEarlyTerminators > 0)
        {
          snprintf(str, 4096,
                   "%s\n%d of your integral curves terminated because they "
                   "reached the maximum number of steps.  This may be indicative of your "
                   "time or distance criteria being too large or of other attributes being "
                   "set incorrectly (example: your step size is too small).  If you are "
                   "confident in your settings and want the particles to advect farther, "
                   "you should increase the maximum number of steps."
                   "  Note that this message does not mean that an error has occurred; it simply "
                   "means that VisIt stopped advecting particles because it reached the maximum "
                   "number of steps. (That said, this case happens most often when other attributes "
                   "are set incorrectly.)\n", str, numEarlyTerminators);
        }
    }

    if (issueWarningForCriticalPoints)
    {
        SumIntAcrossAllProcessors(numCritPts);
        if (numCritPts > 0)
        {
            snprintf(str, 4096, 
                     "%s\n%d of your integral curves circled round and round a critical point (a zero"
                     " velocity location).  Normally, VisIt is able to advect the particle "
                     "to the critical point location and terminate.  However, VisIt was not able "
                     "to do this for these particles due to numerical issues.  In all likelihood, "
                     "additional steps will _not_ help this problem and only cause execution to "
                     "take longer.\n", str, numCritPts);
        }
    }

    if (issueWarningForStepsize)
    {
        SumIntAcrossAllProcessors(numStepSize);
        if (numStepSize > 0)
        {
            snprintf(str, 4096, 
                     "%s\n%d of your integral curves were unable to advect because of the \"stepsize\".  "
                     "Often the step size becomes too small when appraoching a spatial "
                     "or temporal boundary. This especially happens when the step size matches "
                     "the temporal spacing. This condition is referred to as stepsize underflow and "
                     "VisIt stops advecting in this case.\n", str, numStepSize);
        }
    }

    if (issueWarningForStiffness)
    {
        SumIntAcrossAllProcessors(numStiff);
        if (numStiff > 0)
        {
            snprintf(str, 4096, 
                     "%s\n%d of your integral curves were unable to advect because of \"stiffness\".  "
                     "When one component of a velocity field varies quickly and another stays "
                     "relatively constant, then it is not possible to choose step sizes that "
                     "remain within tolerances.  This condition is referred to as stiffness and "
                     "VisIt stops advecting in this case.\n", str, numStiff);
        }
    }

    if( strlen( str ) )
    {
        snprintf(str, 4096, 
                 "\n%s\nIf you want to disable any of these messages, "
                 "you can do so under the Advanced tab.\n", str);

        avtCallback::IssueWarning(str);
    }
}


// ****************************************************************************
//  Method: avtLCSFilter::CreateCacheString
//
//  Purpose:
//      A routine that calculates a string for caching that encodes all the
//      parameters of the LCS: bounds, integration time, and variable name.
//
//  Programmer: Allen Sanderson
//  Creation:   September 5, 2013
//
//  Modifications:
//
// ****************************************************************************

std::string
avtLCSFilter::CreateCacheString(void)
{
  const double* velocitySource = atts.GetVelocitySource();

  std::ostringstream os;

  // Note the string is a brute force string. That is all attributes
  // are munged together regarless of whether they are used. Overly
  // simple but most attributes can not be changed unless an option is
  // selected. As such, for the most part the brute force approach is
  // acceptable.
  os << PAR_Rank() << "  "
     << cycleCached << "  "
     << timeCached << "  "

     << CreateCacheStringSource()

     << atts.GetAuxiliaryGrid() << " "
     << (atts.GetAuxiliaryGrid() ? atts.GetAuxiliaryGridSpacing() : 0) << " "
    
     << atts.GetIntegrationDirection() << " "

     << atts.GetFieldType() << " "
     << atts.GetFieldConstant() << " "
     << velocitySource[0] << " "
     << velocitySource[1] << " "
     << velocitySource[2] << " "

     << atts.GetIntegrationType() << " "
     << atts.GetMaxStepLength() << " "
     << atts.GetLimitMaximumTimestep() << " "
     << atts.GetMaxTimeStep() << " "
     << atts.GetRelTol() << " "
     << atts.GetAbsTolSizeType() << " "
     << atts.GetAbsTolAbsolute() << " "
     << atts.GetAbsTolBBox() << " "

     << atts.GetOperationType() << " "
     << atts.GetEigenComponent() << " "
     << atts.GetEigenWeight() << " "
     << atts.GetOperatorType() << " "
     << atts.GetCauchyGreenTensor() << " "
     << atts.GetClampLogValues() << " "

     << atts.GetTerminationType() << " "

     << atts.GetTerminateByTime() << " "
     << (atts.GetTerminateByTime() ? atts.GetTermTime() : 0) << " "

     << atts.GetTerminateByDistance() << " "
     << (atts.GetTerminateByDistance() ? atts.GetTermDistance() : 0) << " "

     << atts.GetTerminateBySize() << " "
     << (atts.GetTerminateBySize() ? atts.GetTermSize() : 0) << " "

     << atts.GetMaxSteps() << " "

     << atts.GetThresholdLimit() << " "
     << atts.GetRadialLimit() << " "
     << atts.GetBoundaryLimit() << " "
     << atts.GetSeedLimit() << " "

     << atts.GetPathlines() << " "
     << atts.GetPathlinesOverrideStartingTimeFlag() << " "
     << atts.GetPathlinesOverrideStartingTime() << " "
     << atts.GetPathlinesCMFE() << " "

     << atts.GetParallelizationAlgorithmType() << " "
     << atts.GetMaxProcessCount() << " "
     << atts.GetMaxDomainCacheSize() << " "
     << atts.GetWorkGroupSize() << " "

     << atts.GetIssueAdvectionWarnings() << " "
     << atts.GetIssueBoundaryWarnings() << " "
     << atts.GetIssueTerminationWarnings() << " "
     << atts.GetIssueStepsizeWarnings() << " "
     << atts.GetIssueStiffnessWarnings() << " "
     << atts.GetIssueCriticalPointsWarnings() << " "
     << atts.GetCriticalPointThreshold() << " ";

  return os.str();
}

std::string
avtLCSFilter::CreateCacheStringSource(void)
{
  const int*    resolution = atts.GetResolution();
  const double* startPosition = atts.GetStartPosition();
  const double* endPosition = atts.GetEndPosition();

  std::ostringstream os;

  os << atts.GetSourceType() << " ";
  
  if( atts.GetSourceType() == LCSAttributes::NativeMesh )
  {
    os << 0 << " " << 0 << " " << 0 << " "
       << 0 << " " << 0 << " " << 0 << " " << 0 << " "
       << 0 << " " << 0 << " " << 0 << " " << 0 << " ";
  }
  else //if( atts.GetSourceType() == LCSAttributes::RegularGrid )
  {
    os << resolution[0] << " " << resolution[1] << " " << resolution[2] << " ";

    os << atts.GetUseDataSetStart() << " ";

    if( atts.GetUseDataSetStart() ) // Subset start
    {
      os << startPosition[0] << " "
         << startPosition[1] << " "
         << startPosition[2] << " ";
    }
    else // Full
    {
      os << 0 << " " << 0 << " " << 0 << " ";
    }

    os << atts.GetUseDataSetEnd() << " ";
         
    if( atts.GetUseDataSetEnd() ) // Subset end
    {
      os << endPosition[0] << " "
         << endPosition[1] << " "
         << endPosition[2] << " ";
    }
    else // Full
    {
      os << 0 << " " << 0 << " " << 0 << " ";
    }
  }
  
  return os.str();
}

// ****************************************************************************
//  Method: avtLCSFilter::GetCachedDataSet
//
//  Purpose:
//      Checks to see if we have already calculated the LCS.  Works for
//      both resampling and native options (realized through sub-routine calls).
//      This routine uses collective communication to decided whether
//      it can use a cached data set.
//      - For the native resolution, all domains on all MPI tasks must have
//      the LCS cached.
//      - For the resampled version, one MPI task somewhere must find it.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

avtDataTree_p
avtLCSFilter::GetCachedDataSet()
{
    avtDataTree_p rv = NULL;

    if (atts.GetSourceType() == LCSAttributes::NativeMesh)
    {
        int looksOK = 1;

        rv = GetCachedNativeDataSet(GetInputDataTree());
      
        if ((*rv == NULL) && (*(GetInputDataTree()) != NULL))
            looksOK = 0;

        // Gather operation from all procs, if any fails, we all fail.
        looksOK = UnifyMinimumValue(looksOK);

        if (looksOK == 0)
            rv = NULL;
    }
    else //if (atts.GetSourceType() == LCSAttributes::RegularGrid)
    {
        rv = GetCachedResampledDataSet();

        int looksOK = (*rv == NULL ? 0 : 1);

        // Gather operation from all procs, if any fails, we all fail.
        looksOK = UnifyMaximumValue(looksOK);

        if (looksOK == 0)
            rv = NULL;
        else if ((looksOK == 1) && (*rv == NULL))
            rv = new avtDataTree();
    }

    return rv;
}

// ****************************************************************************
//  Method: avtLCSFilter::GetCachedNativeDataSet
//
//  Purpose:
//      Checks the cache to see if we have calculated the LCS on this
//      domain previously.  Checks to make sure parameters match: integration
//      time, etc.
//
//  Programmer: Hank Childs
//  Creation:   December 5, 2011
//
// ****************************************************************************

avtDataTree_p
avtLCSFilter::GetCachedNativeDataSet(avtDataTree_p inDT)
{
    if (*inDT == NULL)
        return NULL;

    int nc = inDT->GetNChildren();

    if (nc < 0 && !inDT->HasData())
    {
        return NULL;
    }

    if (nc == 0)
    {
        //
        // there is only one dataset to process
        //
        int domain = inDT->GetDataRepresentation().GetDomain();
        std::string label = inDT->GetDataRepresentation().GetLabel();
        std::string str = CreateCacheString();

        vtkDataSet *rv = (vtkDataSet *)
          FetchArbitraryVTKObject(SPATIAL_DEPENDENCE | DATA_DEPENDENCE,
                                  outVarName.c_str(), domain, -1, str.c_str());

        if (rv == NULL)
            return NULL;
        else
            return new avtDataTree(rv, domain, label);
    }

    //
    // there is more than one input dataset to process
    // and we need an output datatree for each
    //
    avtDataTree_p *outDT = new avtDataTree_p[nc];
    int cc = 0;
    bool badOne = false;
    for (int j = 0; j < nc; j++)
    {
        if (inDT->ChildIsPresent(j))
        {
            outDT[j] = GetCachedNativeDataSet(inDT->GetChild(j));

            if (*(outDT[j]) == NULL)
            {
                badOne = true;
                break;
            }
            else
            {
                ++cc;
            }
        }
        else
        {
            outDT[j] = NULL;
        }
    }

    avtDataTree_p rv = NULL;
    
    // if we don't have LCS for one domain, then just re-calc whole thing
    if (cc && !badOne)
    {
        rv = new avtDataTree(nc, outDT);
    }

    delete [] outDT;
    return rv;
}


// ****************************************************************************
//  Method: avtLCSFilter::GetCachedResampledDataSet
//
//  Purpose:
//      Checks the cache to see if we have calculated the LCS before.  It
//      also ensures that the parameters of previous calculations are the
//      same: integration time, bounds, variable, etc.
//
//  Programmer: Hari Krishnan
//  Creation:   December 5, 2011
//
// ****************************************************************************

avtDataTree_p
avtLCSFilter::GetCachedResampledDataSet()
{
    std::string str = CreateCacheString();
    vtkRectilinearGrid *rv = (vtkRectilinearGrid *)
      FetchArbitraryVTKObject(SPATIAL_DEPENDENCE | DATA_DEPENDENCE,
                              outVarName.c_str(), -1, -1, str.c_str());

    if( rv )
    {
      int dims[3];
      rv->GetDimensions( dims );
            
      int extent[6];
      rv->GetExtent( extent );

      double bounds[6];
      rv->GetBounds( bounds );
    }
    

    if(rv != NULL)
        return new avtDataTree(rv, 0);
    else
        return NULL;
}
